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Abstract. Black hole-neutron star (BHNS) binaries are important sources of gravitational 
waves for second-generation interferometers, and BHNS mergers are also a proposed engine 
for short, hard gamma-ray bursts. The behavior of both the spacetime (and thus the emitted 
gravitational waves) and the neutron star matter in a BHNS merger depend strongly and 
nonlinearly on the black hole's spin. While there is a significant possibility that astrophysical 
black holes could have spins that are nearly extremal (i.e. near the theoretical maximum), 
to date fully relativistic simulations of BHNS binaries have included black-hole spins only 
up to S/M'^=Q.9, which corresponds to the black hole having approximately half as much 
rotational energy as possible, given the black hole's mass. In this paper, we present a new 
simulation of a BHNS binary with a mass ratio q ~ 3 and black-hole spin S'/A/'^=0.97, the 
highest simulated to date. We find that the black hole's large spin leads to the most massive 
accretion disk and the largest tidal tail outflow of any fully relativistic BHNS simulations 
to date, even exceeding the results implied by extrapolating results from simulations with 
lower black-hole spin. The disk appears to be remarkably stable. We also find that the high 
black-hole spin persists until shortly before the time of merger; afterwards, both merger and 
accretion spin down the black hole. 
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1. Introduction 
1.1. Motivation 

Black hole-neutron star (BHNS) mergers are expected to be a major source for the 
upcoming advanced interferometric gravitational-wave detectors (Advanced LIGO, VIRGO, 
and KAGRA [1, 2, 3]), with event rates in the range ~ 0.2 - 300 BHNS/yr [4]. If the 
mergers leave a massive accretion disk around the black hole, they are also a promising 
setup for short-duration gamma ray bursts (SGRB) [5, 6, 7], which might be observed either 
directly or as an "orphan" SGRB afterglow. If matter is ejected during the merger (either 
from the tidal tail or a disk wind), other signals are conceivable, including radio afterglows 
and/or "kilonovae". (See [8] for a discussion of electromagnetic counterparts.) 

It is to be expected that the parameters of BHNS systems will vary widely; the mass and 
spin of the black hole, in particular, can be substantially different from system to system. 
Population synthesis studies [9, 10] and mass measurements of stellar-mass black holes [11] 
both favor black-hole masses of 8-1 IMq, but they are currently unable to constrain the spin 
of the black hole, because this is set mostly by the spin acquired at the hole's formation [9]. 
This spin at formation will depend on the progenitor of the black hole and the complicated 
core collapse and accretion dynamics of the black hole's birth; existing simulations can 
produce a wide range of spins [12, 13]. 

The black-hole spin can strongly affect the outcome of the BHNS merger. Spin decreases 
the innermost stable circular orbit (ISCO) radius and specific angular momentum for 
prograde circular orbits, and it increases them for retrograde orbits. It is expected that 
aligned black-hole spins will increase the mass of the post-merger accretion disk. Numerical 
simulations confirm this expectation, finding that the disk mass increases dramatically with 
black-hole spin [14, 15]. Indeed, for BHNS systems with black-hole mass > lOM©, only those 
with high dimensionless black-hole spin appear capable of producing significant post-merger 
accretion disks [16, 17, 18]. For a low mass system (Mbh ~ 4:.5Mq), our previous simulations 
found that the mass remaining outside of the black hole after merger increased from 5% of 
the neutron-star mass for a system with an initially nonrotating black hole to 39% for a 
system with a pre-merger black-hole dimensionless spin of x = S/M'^ =0.9 [15]. 

It might appear that this previous study has already effectively covered the range of 
possible black-hole spins and the consequent post-merger behavior. However, x turns out 
to be a poor measure of how close a black hole is to extremality. For example, a black hole 
with X = 0.9 has only about half the rotational energy of an extremal black hole of the 
same mass. More importantly for the accretion phase, the ISCO radius shrinks significantly 
between x = 0-9 s^^id x = 1 (from 2.3Mbh to Mbh [19]). Thus, increasing the spin in 
this range can appreciably change the mass of nuclear matter that is able to avoid prompt 
accretion. It can also increase the efficiency of energy extraction from the accretion disk, 
rj = L/Mc^, where L is the luminosity of the disk, M is the accretion rate, and c is the 
speed of light. Taking rj to be the specific binding energy at the ISCO, as for a thin disk 
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with no ISCO shear, the efficiency for an extremal black hole is higher than that of a % = 0.9 
black hole by a factor of 2.8. Finally, we note that for x = 0.9, the ISCO is still outside 
the ergosphere; the ISCO and, perhaps, significant amounts of disk matter, are only found 
inside the ergoregion for x > 0.94. It is conceivable that new behavior is possible for these 
"extreme" cases of accreting x > 0.94 black holes. 

Another important drawback of our previous study was that the disk was evolved for only 
a short time after its formation. A thick accretion disk of nuclear matter with a mass roughly 
10% of the central black hole could show interesting dynamics. If the disk is sufficiently cool, 
self-gravity could make the disk unstable to spiral density modes or fragmentation. For thin 
disks, stability is determined by the Toomre parameter Qt = kCs/ttGS [20, 21], where Cg is 
the sound speed of the disk, k is the epicychc frequency, S is surface density, and G is the 
gravitational constant. Instabilities in thick self-gravitating disks have been investigated 
numerically in both Newtonian [22, 23, 24, 25, 26] and relativistic [27, 28, 29, 30, 31] 
physics. The instabilities catalogued in these studies are often very sensitive to the angular 
momentum profile of the disk, especially the runaway [32, 33] and the Papaloizou-Pringle 
instabilities [34]. Since the disk is thick and non-Keplerian, its angular momentum proffie 
can only be known by modeling its formation — the BHNS merger, in our case. 

Extreme black-hole spin can also profoundly enhance the electromagnetic signal from 
the post-merger system. The ISCO binding energy is, as we have mentioned, a strongly 
nonlinear function of x- Models of neutrino-dominated accretion flows have also shown that 
a high black-hole spin leads to a higher neutrino luminosity, presumably making more energy 
available for a gamma ray burst [35, 36, 37, 38]. There is also much more rotational kinetic 
energy in the black hole that can be extracted via the Blandford-Znajek process (with the 
luminosity L^z ~ X^) [39]. 

1.2. Summary and overview 

To test the limit of high black-hole spin, high post-merger disk-mass BHNS mergers, 
we numerically model an extreme case. For the case that we consider, the pre-merger 
dimensionless black-hole spin is x = 0.97, aligned with the orbital angular momentum. 
The mass ratio is 3:1, corresponding to a low black- hole mass system. The neutron star has 
a low compaction Mns/-Rns = 0.144, corresponding to a neutron-star radius i?NS ~14km for 
a 1.4M0 (-Rns ~12km for a I.2M0 neutron star). This system is not intended to represent 
a typical black hole-neutron star binary. Rather, it illustrates the dynamics possible at one 
extreme of the astrophysically allowable parameter space. Also, the combination of mass 
ratio 3:1 and compaction 0.144 has been particularly well studied by numerical relativity for 
lower spins, from x = to x = 0.9. With this spin and compaction, we can compare to 
these other results and isolate the effects of extreme spin. 

We begin our evolution 5.5 orbits (~ 21 ms) before merger, when the binary is still 
in quasicircular orbit. We follow the system until ~ 27 ms after merger; this is 13 orbital 
periods of the maximum density region of the post-merger accretion disk, long enough to let 
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us clearly see the settling of the accretion disk to near equilibrium. 

As expected, the extreme spin leads to a very large post-merger accretion disk. Indeed, 
a full 60% of the neutron-star matter is able to avoid falling into the black hole during the 
initial plunge and merger. The vast majority of this matter settles into a massive accretion 
disk. We note that the early disk mass is significantly larger (by about 50%) than what 
was found even for a x = 0.9 system, indicating that the disk mass is very sensitive to the 
black-hole spin when the spin is nearly extremal. 

The matter evolution can be divided into phases, each described in more detail below. 
First, the neutron star is tidally disrupted, leading to an outgoing tidal tail and an ingoing 
accretion stream. Second, the accretion stream circles the black hole and collides with itself. 
This shock rapidly heats the matter and disrupts the accretion stream. Third, the disk settles 
into axisymmetric equilibrium, starting in the inner regions and proceeding outward. This 
requires redistribution of the orbital energy and angular momentum to allow fluid elements 
to settle into circular orbit. In the outer disk, the fluid is nonaxisymmetric and has sub- 
equilibrium average angular momentum. The disk accretes into the black hole slowly on a 
timescale of hundreds of milliseconds. 

The rest of this paper is organized as follows. In Sec. 2, we summarize the numerical 
methods we use to construct and evolve initial data for a black hole and neutron star in a 
circular inspiral, with the black hole having nearly extremal spin. Then, in Sec. 3, we present 
and discuss the results of our simulation, including the emitted gravitational waves and the 
behavior of the black hole's apparent horizon (Sec. 3.1), the behavior of the accretion disk 
that forms following the neutron star's tidal disruption (Sec. 3.2), and the unbound material 
ejected by the merger (Sec. 3.3). We briefly conclude in Sec. 4. 

2. Methods 

To construct initial data for a black hole-neutron star binary with nearly extremal black- 
hole spin, we follow the methods of Ref. [40], which are in part motivated by the methods 
for constructing binary-black- hole initial data with nearly extremal spins in Ref. [41]. We 
find that these methods, which had only been tested for spins up to x = 0.9 so far, can 
be applied without any modification to generate initial configurations with significantly 
higher black-hole spins — although obtaining the same accuracy as for lower spins requires 
the use of a finer numerical grid close to the black hole. Here, we consider a black hole 
with initial spin x = 0.970. Because the initial conditions are not perfectly in equilibrium, 
the spin slightly decreases during the initial relaxation, falling to o/Mbh = 0.967 by time 
t/M = IOOM0 = 1.5 ms after the beginning of the simulation (cf. Fig. 2). We use a 
polytrope 

P = «:pr + rp, (1) 

where P is pressure, p is density, and T is a fluid variable related to the physical temperature. 
As in Ref. [15], we choose F = 2 and choose k, so that the compaction of the star is 
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Figure 1. The trajectory of the apparent-horizon center (thin, black curve) and center of 
mass of the neutron star (thick, red curve) before tidal disruption (i.e., during the first « 20 
ms after the start of the simulation). The orbital eccentricity in this case is e = 0.005±0.001. 

Mns/-Rns = 0.144. Finally, we reduce the orbital eccentricity of the binary by iteratively 
solving for its instantaneous angular velocity and inspiral rate [42]. The last iteration, which 
is used throughout this paper, has eccentricity e = 0.005 ± 0.001. The orbital trajectories 
of the black hole and neutron star before tidal disruption (i.e., during the first ~ 20 ms of 
evolution after the start of the simulation, which corresponds to approximately 5.75 orbits), 
are shown in Fig. 1. 

The coupled evolution of Einstein's equations of general relativity and of the relativistic 
hydrodynamics equations is performed with the SpEC code [43], using our two-grid 
method [44]: Einstein's equations are evolved in the generalized harmonic formalism [45], 
using pseudospectral methods and excision of the black hole interior, while the neutron star 
fluid is evolved on a separate finite-difference grid covering only the regions in which matter 
is present. The general relativistic equations of hydrodynamics are evolved in conservative 
form using a second-order finite volume scheme and high-order shock capturing methods: 
the approximate Riemann problem is solved on cell faces using the Weighted Essentially 
Non-Oscillatory (WEN05) reconstruction algorithm [46, 47] and Harten, Lax, and van Leer 
(HLL) fluxes [48]. Both sets of equations are evolved in time using the third-order Runge- 
Kutta method with adaptive choice of the time step. 

To handle the black hole's high spin during the evolution, we find that we can apply the 
same techniques that have previously been applied to achieve fully relativistic simulations of 
merging black hole-black hole (BHBH) binaries with spins of x = 0.97 [49, 50]. In particular. 
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because we do not impose a boundary condition on the excision surface inside the black hole's 
apparent horizon, for the numerical evolution to be well-posed we must dynamically adjust 
the shape and velocity of the excision surface to prevent incoming characteristic fields [51]. 
During the first ~ 10 ms of the evolution, we permitted the excision surface to fall well inside 
the apparent horizon, provided that there was no danger of incoming characteristic fields, 
but after noticing that this choice had the undesirable side effect of a noticeable increase in 
the measured Einstein constraint violation (cf. Fig. 3), for the remainder of the evolution 
we instead chose to adjust the excision surface so that it tracked the size of the apparent 
horizon. 

We simulated three resolutions through inspiral and tidal disruption, evolving ^ 27 ms 
(low, medium resolution) and ~ 11 ms (high resolution) after merger. During the 
inspiral, these resolutions correspond to approximately 57^, 65^, and 72^ spectral gridpoints, 
respectively, and 100^, 120^, 140^ finite-difference gridpoints, respectively. After the neutron 
star begins to tidally disrupt, we change our domain decomposition, using low, medium, and 
high resolutions with approximately 98^, 111^, and 124'^ spectral gridpoints, respectively, 
and 140^, 160^, 180^ finite-difference gridpoints, respectively. From then on, we use spectral 
adaptive mesh refinement to dynamically adjust the resolution of each spectral subdomain. 
(Note that for all evolutions, since the black-hole spin is aligned with the orbital angular 
momentum, when solving the relativistic hydrodynamic equations we only evolve the region 
above the orbital plane, imposing a symmetry condition on the orbital plane and reducing 
the number of gridpoints in the direction normal to the orbital plane by a factor of 2.) 

We note that we do not find it necessary or advantageous to change to a damped 
harmonic gauge condition, as was done in the high-spin BHBH simulations in Refs. [49, 50]. 
Instead, we use the same gauge conditions used in previous two-grid BHNS evolutions at the 
same mass ratio [15]. 

3. Results 

In this section, we present our main numerical results. First, in Sec. 3.1, we show results 
for the spacetime's evolution. Then, we discuss the properties of the resulting accretion disk 
(Sec. 3.2) and ejecta (Sec. 3.3) after the black hole tidally disrupts the neutron star. 

3.1. Spacetime 

To characterize the behavior of the curved spacetime, we begin by examining the mass 
and spin of the black hole as measured on the apparent horizon. Figure 2 shows, for 
each resolution, the black- hole dimensionless spin x '■= S/M^^ (top panel), Christodoulou 
mass Mch := y^M^7+^S^V~(4^^^lJ (middle panel), and irreducible mass M\„ (lower panel) 
as functions of time. Following the initial relaxation, the masses increase slightly (from 
Mii-r = 3.31 and M^h = 4.20 to Mjrr = 3.33 and M^h = 4.21) while the dimensionless spin 
decreases slightly (from x = 0.970 to x = 0.967). The masses and spins then remain constant 
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Figure 2. The evolution of the black-hole spin x '■= S/M'^yh (top panel), Christodoulou 
mass Afchr (middle panel), and irreducible mass M{„ (bottom panel) as a function of time 
(with times t shown relative to time to. 8, the time when the baryonic mass outside the hole 
has fallen to 80% of its initial value) for a black hole-neutron star merger with an initial 
black hole spin of magnitude x ~ 0.97. 



in time to within our numerical accuracy until approximately 20 ms after the start of the 
simulation, when the neutron star begins to disrupt. During the next ~ 1 ms, the masses 
sharply increase, while the dimensionless spin sharply decreases. The hole's masses then 
continue to slowly increase as its dimensionless spin slowly decreases. 

While all three resolutions exhibit the same qualitative behavior, we were unable to 
demonstrate convergence of the black hole's mass and spin through the end of the simulation. 
Shortly after disruption begins, in the high-resolution simulation the black hole experiences 
an anomalously large decrease in the dimensionless spin; through experimentation, we found 
this decrease to be sensitive to the details of dynamic regridding (i.e., the algorithm used to 
resize the finite-difference grid as the region containing matter expands — cf. Sec. II B 1 of 
Ref. [15]). This could be a consequence of the extremely large amount of mass remaining 
in the disk. Apparently, there are certain special distributions of matter on our grid that, 
upon a regrid, result in unexpectedly large interpolation error near the excision surface and 
a corresponding small but abrupt change in the evolution. Since our primary concern in this 
paper is the behavior of the disk and ejecta, we have chosen to leave our efforts to resolve 
the black-hole spin after disruption (perhaps by improving our regridding method) for future 
work. 

Figure 3 further demonstrates our difficulty in obtaining convergence after tidal 
disruption. Before disruption (i.e. during the first 20 ms of the simulation), the normalized 
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Figure 3. The normalized Einstein constraint violation (as defined in Eq. (71) of Ref. [45]) 
as a function of time (with times t shown relative to time to.Si the time when the baryonic 
mass outside the hole has fallen to 80% of its initial value) for low (dashed, black), medium 
(solid, red), and high (dotted, green) resolutions. In each simulation, the initial black- hole 
spin is aligned with the orbital angular momentum and has a magnitude of x = 0.97. The 
constraints rise around time t—to.s = —14 ms because the apparent horizon is moving farther 
from the excision surface. After modifying the simulation's algorithm so that the excision 
surface and apparent horizon would nearly coincide (at approximately t — to. 8 = ^10 ms), 
the constraint violation drops until the neutron star begins to disrupt. 



constraint violation appears to be convergent, as expected. The constraints do rise around 
t — to,s ~ — 14 ms as the excision surface falls far inside the apparent horizon, but this trend 
reverses after a modification of our method to keep the excision surface close to the apparent 
horizon. After interpolating to our higher-resolution domain (around t — to.8 ~ ~2 ms), the 
constraints sharply fall at first, but then sharply spike during the merger. The subsequent 
bump in the constraints corresponds to the small bump in the dimensionless spin and the 
small drop in irreducible mass seen around t — to.s ~ 7 nis in Fig. 2. 

We conclude this section by examining the dominant mode of the gravitational waveform 
(Fig. 4). Because we plot = {cP/dt"^) {h^ — iky) instead of the wave amplitude h, the 
spurious "junk" gravitational radiation emitted during the initial relaxation is clearly visible. 
Each resolution shows a small, secondary burst of gravitational waves after the primary wave 
has terminated, but this burst is not resolved. Further simulations at higher accuracy would 
be needed to verify whether this burst is a physical result of the disruption or whether it is 
a numerical artifact. 
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Figure 4. The real part of the dominant i ^ m = 2 mode of the gravitational waveform 
rMsun^^, where r is the radius of the sphere on which the gravitational waveform is extracted 
(here the outermost sphere, with a radius of r = I54OM0), Mgun = ^Mq, and ^'^ is the 
Newman-Penrose scalar corresponding to outgoing gravitational waves. The inset zooms in 
on the part of the waveform enclosed by the dashed box. 

3.2. Accretion disk 

After tidal disruption, the nuclear matter forms an outgoing tidal tail and an ingoing 
accretion stream. At a time t = to.s + 10 ms (10 ms after the time when 80% of the 
initial baryonic mass remains outside the hole), ~ 23% of the initial baryonic mass has 
moved outward and left the finite-difference grid, and ~ 31% of the initial baryonic mass 
remains on the finite-difference grid (Fig. 5). Of the mass that has left, ~ 17% appears to 
be unbound and quickly leaves the grid (see Sec. 3.3), while another !^ 6% is weakly bound 
but extends far enough from the black hole to leave the finite-difference grid. This material 
would have fallen back onto the disk at later times (~ 10^ ms) and its loss is one reason we 
limit our simulations to the early to middle (< 30 ms) post-merger evolution. 

During this tidal disruption phase, the matter remains cold and degenerate. The 
infalling stream of matter circles around the black hole and collides with itself, producing 
an intense shock that travels through and disrupts the incipient accretion flow. Nearly all of 
the entropy generation in the evolution occurs during this crucial millisecond, which could 
be labeled the shock phase of the evolution. Shock heating renders the matter completely 
nondegenerate, and hot matter flows to the outer disk region (r > 100km) over the next 
~ 10ms. By the end of this time, a clear hierarchy of velocities establishes itself through 
the disk out to 200km. (See Fig. 8.) Inside this region, the azimuthal velocity dominates 
and is roughly twice the sound speed throughout (fir ^ 2cs). The orbital period sets the 
dynamical timescale; we may define a settling radius rscttic(^) as the radius at which one 
orbital period has passed since merger. By the end of our evolution (30 ms postmerger). 
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Figure 5. The evolution of the baryonic mass Mb as a function of time t for the 
medium resolution black hole-neutron star merger with an initial aligned black-hole spin 
of magnitude x — 0.97. The solid, green line is the total baryonic mass Mb outside the 
black hole, the dashed, black line is the total rest mass Afgiid remaining on the finite- 
difference computational grid, and their difference Mout is the dotted, red line. All masses 
are normalized by the baryonic mass Mh_Q at the start of the simulation. Times are shown 
relative to time to. 8; the time when the baryonic mass outside the hole has fallen to 80% of 
its initial value. 



'^settle = 190 km. Outside rsettie, the disk does not have time to reach equihbrium. The 
radial velocity is nonnegligible in the outer disk, but the radius at which this infall becomes 
supersonic recedes with rsettie- 

In Figure 6, we show the surface density S at three times during the settling phase. 
The density peaks at a circumferential radius of around 30km; the location of the peak 
remains fixed, while the overall density slowly decreases everywhere as the mass is depleted 
by accretion. Also shown in this plot is the RMS deviation in S at each radius due to 
deviations from axisymmetry. We see that the inner disk is quite axisymmetric, but that 
deviations from axisymmetry reach order unity at around rsettie(^)- A Fourier analysis of S 
at late times shows that the deviation from axisymmetry primarily subsists in trailing spiral 
features, although localized twists in which a mode is leading over a small radial range are 
occasionally seen in some modes. Inside 100km, 5S/Sm=o ~ 10~^, but the nonaxisymmetry 
increases with radius and reaches order unity at around rgettic; with m = 1 dominating in 
this region. Only in the inner region can these modes be regarded as linear perturbations 
of an equilibrium system; here they provide a natural way to understand the extraction of 
angular momentum from the inner disk which must occur if the low- angular-momentum 
outer disk (see below) is to settle. In addition to these smooth modes, equilibrium in 
'^settie/S < r < Tsettie IS sometimes disturbed by weak, localized sharp features, presumably 
weak shocks where gas on eccentric orbits hits the more settled disk. These features are 
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Figure 6. The density profile of the accretion disk, shown at times t =8, 17, and 26ms after 
merger. For each time, two lines are shown. The top line is the vertically integrated surface 
density S at each radius, averaged over all angles. The bottom line is the rms deviation 
from the average due to nonaxisymmetry. A dot on the bottom line indicates the location 
of Tscttic at this time. In this and all other radial plots, the radius shown is constructed from 
the proper circumference at the equator. 

most clearly seen in the pressure force (see below). 

In Figure 7, we plot the specific orbital energy {E = Ut) and angular momentum 
(L = —u^/ut) on the equator as a function of radius. For comparison, we plot the expected 
curves for circular orbit geodesies in this spacetime metric. The geodesic E and L are the 
expected orbital parameters for a stationary disk taking into account the disk's self-gravity 
(which is included in the numerical metric) but not pressure forces. Except in the very inner 
disk, we see significant deviations from geodesic orbits, as would be expected for a thick disk. 
Figure 7 also includes E and L for equilibrium circular orbits, i.e. stationary circular orbits 
taking into account the pressure force for the given fiuid profile. These agree with the actual 
E and L up to about 100km. Beyond this, the energy curves continue to agree (accounting 
for deviations due to localized shocks), but the angular momentum profiles fiatten more 
rapidly than the equilibrium curve (although the latter seems to become somewhat noisy) 
and become sub-equilibrium. The gas at these radii was ejected (from the tidal tail or the 
post-merger shock) into eccentric orbits, as seen from the L deficit and the nonzero radial 
velocity (cf. Fig. 8). Equivalently, the nonaxisymmetric modes carry net negative angular 
momentum. Pressure support becomes especially important in the outer disk, as can be seen 
from the deviation between geodesic and equilibrium curves. We also have observed that the 
disk thickness, height H divided by radius r, increases with radius from about H/r 0.2 in 
the inner disk to about H/r ^ 0.35 in the outer disk. 

There are no clear signs of instability or turbulence in the disk, a fact that itself should 
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be explained. The lack of obvious global corotation instability of the kind found in some 
other massive disk simulations is a natural consequence of the angular momentum profile. 
Instability is only expected for L oc r", < 2 — [34] with growth rate decreasing rapidly 
as n increases above zero, and our disk has n ~ 0.3 in the equilibrated region. This region is 
also stable against the effects of self-gravity, with Toomre parameter no lower than Qt ~ 400. 
This stability is a consequence both of the shear and the thermal pressure. Radiative cooling 
may later remove some of the disk's heat, but even removing all of it would leave Qt ~ 10. 
Shear and thermal energy also conspire to protect the disk from convective instability. In 
Newtonian physics, the Solberg-Hoiland stability criteria are 



where p is the density, S the entropy, P the pressure, and Cp the heat capacity. The 
relativistic version of these equations is slightly more complicated but basically similar [52]. 
Figure 9 shows the average specific entropy as a function of radius during the settling phase. 
The entropy has a minimum close to the density maximum, while the gas on the edges is 
hotter. Thus VP ■ \/S < except near a small ring near r = 40km, and this feature is 
stabilized by the strong shear. In fact, the epicyclic frequency is much greater than the 
Brunt- Vaisala frequency out to r ^ 80km. Beyond this, the S gradient steepens, the L 
gradient flattens, and buoyancy becomes an important restoring force. Vertical convection is 
also suppressed by a strong positive VzS. The disk will be unstable to the magnetorotational 
instability, since dQ"^ /dlnr + N'^ < throughout, but this does not appear in our simulations, 
which do not include magnetic fields. 

The disk evolves slowly due to the settling of the outer disk and accretion of gas in 
the inner disk (the two perhaps connected by transfer of angular momentum). The late- 
time accretion rate is M 2Mq s~^, giving a lifetime of r = M/M ~ 200ms. Accretion 
increases the spin and mass of the black hole, but since the accreted matter has low angular 
momentum [L ^ 10, see Fig. 7), the total effect on the dimensionless spin x = S/M"^ will 
be small even for accretion of most of the disk mass. (At late times in our evolution, x is 
slowly decreasing.) 

Our simulation does neglect several effects that could influence the accretion disk's later 
evolution. A magnetic field could produce a jet, powered either by the accretion flow or 
the high black hole spin. It would certainly render our disk magnetorotationally unstable. 
Although we find that MRI turbulence is not needed for strong early-time accretion, it will 
still alter the disk evolution in important ways. Second, neutrino radiation will cool the gas 
and make it less neutron-rich. Assuming a simple ideal gas plus radiation pressure equation 
of state, we infer that the temperature in the disk is of order lOMeV. (In this model, the 
disk is gas pressure-dominated everywhere except the inner edge.) Using a simple diffusion 
approximation, one can estimate a luminosity L^, ~ lO^^erg s~^, which could deplete the 
disk's thermal energy on a timescale of ~ lO^ms. More likely, the cooling will balance 
the dissipative heating introduced by magnetic turbulence. Neutrino radiation could itself 
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Figure 7. The specific orbital energy E and specific angular momentum L at three times 
after merger. For tlie final time, we include two other curves for comparison. First, we 
include E and L for circular orbit geodesies in the given metric (which will include the 
disk's self gravity). Second, we show E and L for equilibrium circular orbits given the 
current metric and the pressure forces from the current pressure and density profiles. For 
the equilibrium energy, we take the median in angle rather than the mean. The mean would 
reveal spikes at 140km and 190km due to local shock features. The geodesic curves for the 
earlier times are nearly identical to those at late times, while the equilibrium curves differ 
only at large radii. 





Figure 8. A break-up of the components of the fluid velocity, shown at three times after 
merger, as a fraction of the speed of light. Also shown is the sound speed Cg- The extents on 
the horizontal and vertical axes are the same in all three panels. The radial velocity itself 
is negative (infall). 
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Figure 9. The entropy profile of the accretion disk, shown at times t ~8, 17, and 26ms 
after merger. The entropy for our F-law equation of state is S = log( p(p'^) )■ For each 
radius, a density-weighted average is carried out over the and z coordinates. 



induce unstable entropy or composition gradients and drive convection. Finally, our grid 
limitations have forced us to allow gravitationally bounded, outward-moving material to leave 
the finite-difference grid. This matter should eventually fall back onto the disk, perturbing 
and perhaps shocking it. Tracking the fallback matter will be a computational challenge for 
future simulations that wish to evolve disks like this one further in time than was done here. 

We should note that the amount of matter remaining outside of the black hole after 
merger is not only larger than in any previous studies of BHNS mergers: it is also significantly 
above the predictions obtained by extrapolating results for lower spin black holes in the high 
spin regime. The remnant mass for lower spin systems is generally well approximated by [18] 



with -Risco the radius of the innermost stable circular orbit around the black hole. Figure 10 
shows the predicted remnant mass for BHNS mergers at g = 3, Mns/-Rns = 0.144, together 
with numerical results at x = (0,0.5,0.75,0.9,0.97) [14, 17]. We see that the analytical 
predictions perform well for low-spin black holes (x < 0.9), while strong deviations are 
visible for the higher spin configurations. This is not surprising: the analytical results are 
only supposed to be valid for relatively low mass remnants (Mrem < 0.2Mq). But it does 
emphasize the need for further studies of BHNS mergers at high black-hole spin in order to 
properly model the characteristics of the post-merger remnant in that regime — especially 
considering that the most energetic SGRBs observed might require the presence of such 
high-mass accretion disks [53]. 
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Figure 10. The evolution of the baryonic mass for black hole-neutron star mergers with 
different initial black-hole spins of magnitude x aligned with the orbital angular momentum. 
The mass for x = 0.97 is shown for the three different resolutions (low, medium, high) 
simulated. The masses from mergers with X = 0; X = 0-5, and x = 0.9 are taken from Fig. 5 
of Ref. [15]. Top panel: The baryonic masses are shown as a function of time t — tg.s, where 
to. 8 is the time at which the baryon mass outside the hole has fallen to 80% of its initial value. 
Bottom panel: The baryonic mass at time t = tp.s + 10 ms (black crosses) as a function of 
the initial black-hole spin. The medium resolution of the x = 0.97 simulation is shown. 
Relative errors of 15% are used for all simulations, an estimate obtained by extrapolating 
our results to infinite resolution while assuming second order convergence of the code. For 
comparison, also shown are i) the final disk mass obtained by the UIUC group in Fig. 13 of 
Ref. [14] for case "B" and ii) the theoretical prediction derived in Ref. [18]. 



3.3. Ejecta 

In addition to the formation of an accretion disk, the tidal disruption of the neutron star 
can lead to the ejection of a significant amount of unbound material. This material is 
initially neutron rich, and its radioactive decay can lead to detectable optical afterglows 
("kilonovae" [54, 8]), as well as contribute to the formation of heavy elements (r-process 
nucleosynthesis). The deceleration of the ejecta in the interstellar medium can also cause 
emission in the radio band [8]. For binary neutron stars [55] and low-spin BHNS [56] binaries, 
general relativistic simulations show that only a small fraction of the neutron star material 
(< 1%) is unbound, with kinetic energies of 10^^ ergs — 10^^ ergs [55]. 

Measurements of the unbound mass are fairly inaccurate in our code, especially for very 
energetic and massive ejecta. This is due to the fact that the fluid equations are only evolved 
up to 20Mbh — 30Mbh from the center of mass of the system, while measurements of the 
properties of the ejecta are best done far away from the black hole and in low-density tidal 
tails (the condition < —1, used to determine whether material is unbound, is only valid 
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for a stationary metric and pressureless fluid). Additionally, the accuracy of the evolution 
far away from the black hole is lower than in the forming disk, where most grid points 
are concentrated. The development of techniques allowing us to efficiently evolve the fluid 
equations at larger distances without losing accuracy close to the black hole (e.g. adaptive 
mesh refinement) will be required to avoid these issues. In a recent paper [57], we showed 
that the main source of error when measuring the mass of the ejecta is generally the large 
grid spacing in the far zone. The same holds in this simulation: we find 

Mej = 0.26Mq ± O.I6M0 (5) 

assuming 2nd order convergence between the medium and high resolution. This confirms 
that BHNS mergers with rapidly rotating black holes are significantly more favorable to the 
ejection of neutron-rich material. In [57], we found that for higher mass ratio systems [q = 7) 
and lower spins x = 0.9, a few percent of a solar mass is hkely to be unbound, already a 
noticeable improvement compared to low-spin systems. The ejecta in the case studied here, 
which is of course an extremely favorable configuration, is an order of magnitude larger — 
thus showing that for rapidly spinning systems, it is conceivable that a large fraction of the 
mass remaining outside of the black hole at late time is unbound. 

Measurements of the kinetic energy of the ejecta are even less reliable than its mass. 
Indeed, the kinetic energy is dominated by the most relativistic parts of the ejecta, which are 
also the most poorly resolved. From measurements of Ut as the outward-moving material is 
leaving the finite-difference grid, we estimate that the ejected material has a median velocity 
V ~ 0.5c, and kinetic energy E'ej ~ 10^^ ergs — 10^^ ergs. These numbers should, however, 
be considered only as order of magnitude estimates. Nonetheless, this is an extremely large 
amount of energy, which would cause the emission of a detectable radio signal for a large 
fraction of the mergers within the range of Advanced LIGO as the ejecta slows down in the 
interstellar medium. 

4. Conclusion 

We have simulated the merger of a black hole-neutron star system with a premerger black- 
hole spin parameter of x =0.97, the highest spin yet attempted for modeling such systems. 
Given the strongly nonlinear dependence of many aspects of the merger on x (black-hole 
rotation energy, ISCO location, accretion efficiency), this simulation carries the numerical 
exploration of BHNS mergers into a distinctly new regime not sampled even by our previous 
simulations with x =0.9. These nonlinearities manifest themselves clearly in our results. 

Upon disruption of the neutron star, less than half of the nuclear matter is promptly 
accreted into the black hole. Almost 20% appears to be ejected from the system in an 
unbounded tidal tail outflow, and the rest settles into a massive accretion disk. Both the 
ejecta and disk masses are higher than any previous fully relativistic BHNS simulation. 
This in itself is unsurprising; both theoretical considerations and previous simulations would 
lead us to expect that in higher-x systems more matter should evade prompt accretion. 



Disk formation in the disruption of a neutron star by a nearly extremal black hole 17 

What is notable is that the disk mass exceeds even the expectations for this spin based on 
extrapolating the trends of lower-spin systems. If such systems occur in nature, they must 
be particularly spectacular events even by the standard of BHNS mergers. 

We have closely followed the behavior of the horizon and of the accretion disk. Both 
merger and accretion have the net effect of decreasing the dimensionless spin x of the black 
hole. While the disk is thick and self-gravitating, it appears to be quite stable, settling 
to an axisymmetric quasistationary state and evolving only slowly under the influences of 
outer-disk settling and accretion-induced mass decrease. Its subsequent evolution will be 
driven in part by physical processes not included in these simulations. 

That existing numerical relativity techniques can successfully treat such an extreme 
system essentially unaltered is encouraging. Further studies of the extreme-spin regime 
of BHNS parameter space can now be attempted; particularly interesting will be the 
exploration of systems with higher black-hole mass. These are expected to be more common 
astrophysically, and they may also be more intriguing from a relativist's point of view, 
since the disruption event happens nearer to the horizon, as measured by the disruption 
separation divided by Mbh- The investigation of the accretion system studied in this paper 
should also be completed by incorporating the remaining crucial physics, particularly the 
magnetic turbulence. 
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